Calcolo numerico avanzato e Grafica¶
Non re-inventare la ruota!¶
Python (come anche Matlab) mette già a disposizione numerose librerie (Toolbox per Matlab) per il calcolo scientifico e la rappresentazione grafica dei risultati.
- NumPy : calcolo scientifico, algebra, statistica,...
- SciPy : calcolo scientifico, algebra sparsa, integrazione di equazioni differenziali, minimizzatori,...
- Scikit-learn : Machine Learning e Pattern Recognition
- Scikit-image : Image Processing
- Matplotlib : plot e grafica
In Matlab tutte queste funzioni le ritroviamo già nei Toolbox di default (ad eccezione dell'image processing che necessita di un Toolbox a parte...facilmente scaricabile)!
Vedere i propri dati¶
Il primo step da eseguire quando si hanno dei dati è riuscire a visualizzarli nel modo opportuno.
Per farlo Python ci mette a disposizione numerose funzioni, per permetterci di rappresentare i nostri dati in diversi modi ed in modo estremamente facile!
Supponiamo di avere un set di punti sperimentali in più dimensioni...
import numpy as np
data = np.random.randn(1000, 2) # matrice 1000x2 di numeri random estratti da una distribuzione normale
Avendo dati multi-dimensionali per prima cosa vediamo di graficarli lungo ogni dimensione...
import matplotlib.pylab as plt
fig, (ax1, ax2) = plt.subplots(1, 2) # dichiaro una figura(fig) con due assi di plot(ax1, ax2)
ax1.plot(data[:, 0], marker = 'o', markersize=3, color = 'k', label='1st dimension') # plot della prima colonna di dati
ax1.set_xlabel("Index", fontsize=14) # label delle ascisse
ax1.set_ylabel("Data", fontsize=14) # label delle ordinate
ax1.legend(loc='best', fontsize=14) # legenda
ax2.scatter(np.arange(0, len(data)), data[:,1], marker = '*', color='r', label='2nd dimension') # scatter plot
ax2.set_xlabel("Index", fontsize=14)
ax2.set_ylabel("Data", fontsize=14)
ax2.legend(loc=3)
plt.show()
Ma io avevo generato dati gaussiani?!¶
Per vedere se effettivamente sono "guassiani" dovrò crearne un istogramma!
plt.figure() # creazione di una figura
plt.hist(data[:,0], color='k', label='1st dim', alpha = .5) # creazione di un istogramma
plt.hist(data[:,1], color='r', label='2nd dim', alpha = .5)
plt.xlabel("x")
plt.ylabel("y")
plt.legend(loc='best')
plt.show()
Posso però anche guardare le due dimensioni insieme con un bel plot 2D
plt.figure()
plt.hist2d(data[:,0], data[:,1], bins = 30) # creazione di un istogramma 2D
plt.xlabel("1st dim")
plt.ylabel("2nd dim")
plt.colorbar() # creazione della colorbar
plt.show()
Vogliamo esagerare?! Facciamolo anche in 3D!¶
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure()
ax = fig.gca(projection='3d')
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y) # creazione della griglia di valori per le coordinate x,y
mean = np.array([0., 1.]) # vettore di media
Cov = np.array([[ 1. , -0.5], [-0.5, 1.5]]) # matrice di covarianza
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y
# This einsum call calculates (x-mu)T.Sigma-1.(x-mu) in a vectorized way across all the input variables.
fac = np.einsum('...k,kl,...l->...', pos-mean, np.linalg.inv(Cov), pos-mean)
Z = np.exp(-fac / 2) * np.sqrt((2*np.pi)**len(mean) * np.linalg.det(Cov))
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False) # creazione del plot 3D
plt.show()
Ora un po' di conti però..¶
La lezione scorsa abbiamo iniziato a vedere alcune delle funzionalità di NumPy, legate all'istanziazione di oggetti ed alle funzioni base per il calcolo vettoriale.
SciPy integra molte delle funzionalità di NumPy e le arricchisce con ulteriori algoritmi avanzati. Iniziamo a vederne qualcuno...
Algebra lineare - LinAlg¶
All'interno di NumPy troviamo il modulo numpy.linalg con diversi utili algoritmi per l'algebra lineare. Allo stesso modo scipy.linalg ne fornisce altri. Tuttavia alcune funzioni di base le troviamo a parte!
Iniziamo con qualcosa di facile...l'algebra tra vettori e matrici!
Abbiamo visto che scrivere
c = a * b
restituisce il prodotto dei singoli elementi (Hadamard product). E allora il prodotto scalare?!
import numpy as np
import scipy as sp
a = np.array([[1,2,3,4,5]])
b = np.array([1,2,3,4,5])
c = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15]])
print (np.dot(b, b)) # prodotto scalare tra b e b
print (sp.dot(b, b))
print (np.dot(c, a.T)) # prodotto scalare tra c e a-trasposto
print (np.dot(a, c.T))
print (np.inner(b, c)) # prodotto scalare tra b e c
print (np.inner(c, b))
print (c.transpose() == c.T)
ATTENZIONE!¶
Mentre in Matlab ogni array è automaticamente considerato un vettore in senso matematico, in Python un semplice array ha solo 1 dimensione! Per questo la funzione dot è un po' più noiosa! Per evitare problemi dovuti alla trasposizione si può utilizzare la funzione inner per fare appunto l'inner product tra due vettori. Occhio però a cosa volete fare!!!
Per tutto il resto c'è la trasposizione .transpose() o (per gli amici) .T.
NOTA: se vi sentite più forti in matematica che in programmazione, senza pensare ai vari problemi di trasposizione potete utilizzare sempre la funzione np.einsum!
Vediamo adesso qualcosa di un po' più avanzato e "difficile". Il calcolo degli autovalori ed autovettori.
Remember, remember...¶
A livello pratico non sempre vi serviranno tutti gli autovalori/autovettori nei vostri calcoli. Molto spesso servono solo i maggiori o minori. Ottimizzazione e scrittura smart di un codice vuol dire anche evitare gli eccessi.
Per questo bisogna fare un po' di attenzione! Soprattutto se lavorate con matrici con tanti valori nulli!
Altra cosa da tenere sempre a mente è l'ordinamento degli autovalori/autovettori!
a = np.eye(5)+np.random.rand(5,5) # matrice data dall'identità sommata ad un rumore uniforme
eigenval, eigenvec = np.linalg.eig(a) # funzione per autovalori e autovettori di una matrice
print ("autovalori = ", eigenval)
print ("autovalori ordinati = ", np.sort(eigenval)[::-1]) # sort degli autovalori in ordine decrescente
eigenval, eigenvec = sp.sparse.linalg.eigs(a, k = 2) # calcolo dei primi 2 autovalori
print ("primi due autovalori = ", eigenval)
print ("primi due autovalori ordinati = ", np.sort(eigenval)[::-1])
A chi piace fare i conti a mano...¶
- Vi piace fare i conti a mano?
- Avete in mente una bella formula per trovare il risultato che cercate e non avete voglia di scrivervi un algoritmo?!
- Vi basta vedere la soluzione in formula per essere contenti?
Nessun problema!: Python e Matlab supportano anche il calcolo simbolico!
import sympy
x = sympy.Symbol('x') # dichiarazione del simbolo 'x'
y = 1 + x + x**2 # funzione di simboli
subs_dict = {x : 2} # sostituisco 2 ad x
y.subs(subs_dict) # ogni volta che incontro x sostituisco 2
print (sympy.solve(y, x, dict=True)) # restituisce in un dizionario le possibili soluzioni di y con x = 2
y.diff(x) # derivata di y rispetto ad x
Equazioni differenziali¶
Un altro importante argomento per un fisico è l'integrazione di equazioni differenziali.
Ci sono tantissimi metodi più o meno accurati per risolverle, utilizzando schemi di integrazioni ad accuratezza crescente (Eulero, RK2, RK4, ...)
Proprio per questo motivo SciPy mette a disposizione un intero modulo scipy.integrate.
L'equivalente di Matlab sono le varie funzioni ode da scegliere accuratamente in relazione alla stabilità del problema.
from scipy.integrate import odeint
Facciamo un esempio:
$ \dot{y} = -\alpha y $
Tutti noi sappiamo che questa equazione, data la condizione iniziale $y(0) = y_0$, avrà per soluzione un esponenziale in funzione di $\alpha$.
Non ci resta che integrarla e vedere con i nostri occhi che questo sia giusto!
def derivata(y, t):
return -y
time = np.linspace(0, 10, 20)
y0 = 10.0
yt = odeint(derivata, y0, time) # integrazione dell'equazione differenziale: (func, condizione iniziale, ascisse)
import matplotlib.pylab as plt
plt.figure()
plt.scatter(time, yt, marker = 'o', color = 'k', label = 'exponential')
plt.xlabel('time', fontsize = 14)
plt.ylabel('y(t)', fontsize = 14)
plt.title('Exponential ODEint', fontsize = 16)
plt.legend(loc = 'best')
plt.show()
Supponiamo di voler però verificare che la soluzione sia proprio un esponenziale.
Per prima cosa dovremo avere dei punti sperimentali un po' più "realistici".
Poi dovremo fittare la nostra funzione.
Mettiamoci all'opera allora..
y_obs = yt.ravel() + np.random.rand(len(yt))*2-1 # dati teorici sommati ad un rumore uniforme tra [-1,1]
plt.plot(time, y_obs, 'bo', label = 'sperimentali')
plt.plot(time, yt, 'k-', label = 'teorici')
plt.legend(loc = 'best')
plt.show()
from scipy.optimize import curve_fit
def decadimento(tempo, alfa):
return 10 * np.exp(-alfa*tempo)
fit_alfa, var_alfa = curve_fit(decadimento, time, y_obs, p0=[0.9]) # (func, ascissa, ordinata, condizione iniziale parametri)
std_alfa = np.sqrt(var_alfa)
plt.plot(time, yt, color='r', linewidth=8, alpha = 0.3)
plt.plot(time, y_obs, 'o', label='dati sperimentali');
y_hat = decadimento(time, fit_alfa)
plt.plot(time, y_hat, '-r', label='curva fittata');
plt.legend(loc = 'best', numpoints = 1)
plt.show()
Image Processing¶
Per l'Image Processing in Python una libreria abbastanza completa è data da Scikit-image, nella quale troviamo le principali funzioni di sogliatura e filtering.
Per il plotting delle immagini continueremo ad utilizzare la nostra Matlplotlib.
Remember, remember:¶
Un'immagine è pur sempre una matrice per cui tutto quello visto fin'ora continua a valere anche per l'image processing!!
Quando si testano nuovi algoritmi di Image Processing o Machine Learning è utile avere dei dati a disposizione. Ovviamente in rete potete trovare tutti i dataset che volete.
Senza fare tanta strada però, installando Python avete anche "inconsciamente" scaricato quelli che sono i data sets "standard".
Vediamo intanto quello per le immagini...
from skimage import data, io, filters
image = data.coins()
plt.figure()
plt.imshow(image) # funzione per il plot di immagini/matrici
plt.show()
print (image.shape)
print (image.dtype)
print (image.max())
print (image.min())
E se adesso volessi fare un po' di processing?!¶
Ad esempio potremmo individuare gli edges nell'immagine
edges = filters.sobel(image)
plt.figure()
io.imshow(edges)
io.show()
E vedere come si comportano vari filtri al variare dei parametri...
from scipy import ndimage as ndi
from skimage import feature
edges1 = feature.canny(image)
edges2 = feature.canny(image, sigma=3)
# display results
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3),sharex=True, sharey=True)
ax1.imshow(image, cmap=plt.cm.gray) # plot con una colormap in scala di grigi
ax1.axis('off')
ax2.imshow(edges1, cmap=plt.cm.gray)
ax2.axis('off')
ax2.set_title('Canny filter, $\sigma=1$', fontsize=20)
ax3.imshow(edges2, cmap=plt.cm.gray)
ax3.axis('off')
ax3.set_title('Canny filter, $\sigma=3$', fontsize=20)
fig.tight_layout()
plt.show()
Machine Learning¶
Finiamo con un po' di Machine Learning!
Uno dei punti fondamentali del Machine Learning sono sicuramente i Classificatori! Ce ne sono tanti, ce ne sono di più precisi e di meno precisi, di più belli e di più brutti...
...ma la cosa più bella di tutte è che ci sono quasi tutti già implementati in Python e Matlab!!!
Per Python faremo riferimento alla libreria Scikit-learn.
Vediamo allora di implementare la ben nota Bayesian Quadratic Discriminant Analysis!
from sklearn import datasets
iris = datasets.load_iris()
data = iris.data
label = iris.target
plt.figure()
plt.scatter(data[:50, 0], data[:50, 1], marker='o', color='k')
plt.scatter(data[50:100, 0], data[50:100, 1], marker = '*', color='r')
plt.scatter(data[100:, 0], data[100:, 1], marker = 'd', color='g')
plt.xlabel("1st dim", fontsize=14)
plt.ylabel("2nd dim", fontsize=14)
plt.show()
from sklearn.qda import QDA
from sklearn.cross_validation import LeaveOneOut
Loo = LeaveOneOut(len(label)) # creazione di un oggetto di iterazione leave-one-out; ritorna la coppia di indici
classifier = QDA() # dichiaro il classificatore
Classes = np.empty(len(label)) # vettore in cui salverò le predizioni del LOO
for train, test in Loo: # ciclo sulle coppie di LOO
Training = data[train, :] # estrazione del data set di training
Test = data[test, :] # estrazione del data set di test
classifier.fit(Training, label[train]) # fitto il training
Classes[test] = classifier.predict(Test) # predico il test
print ("Score QDA %d over %d = %f%%"%(len(Classes[Classes == label]), #numero di predizioni la cui label è uguale a quella vera
len(label),
float(len(Classes[Classes == label]))/len(label)))
ORA TOCCA A VOI !¶
L'obiettivo è quello di riscrivere il classificatore QDA step-by-step! E' vero che non bisogna riscrivere quello che già esiste...però il QDA di Sklearn ogni tanto dà qualche problema! Infatti non sempre riesce ad invertire la matrice di covarianza.
Per risolverlo basta aggiungere un $\epsilon$ agli elementi sulla diagonale della covarianza in modo che la matrice non sia più singolare!
Per chi non se la ricordasse: la formula da implementare è
$$ f(\mathbf{x}) = -0.5 * (\mathbf{x}-\mathbf{\mu_i})^T*\Sigma_i^{-1}*(\mathbf{x}-\mathbf{\mu_i}) -0.5*\log(det(\Sigma_i)) $$dove la $i$ considera le varie classi!
Qualche consiglio...¶
- Con np.unique(label) avete un set delle classi
- Ogni classe deve avere la sua covarianza...quindi ci vorrà un bel ciclo
- la X identifica il TEST
- la classe predetta è scelta in base al valore MAX di f(X)
...sono troppo buono!